T-shaped dimer of coronene 
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An evidence of importance of the T-shaped configuration of coronene dimer is pre- 
sented. That is, the dimer's lowest energy configuration is not necessarily a stack, 
as it might had been expected a priori. This is a surprising result for dimer of 
such a large polycyclic aromatic hydrocarbon (PAH) as coronene (C24H12). The 
energy of the T-shaped configuration at all considered levels of density functional 
theory (B3LYP,PBE/6-31-hG(d),D95,cc-pVDZ,cc-pVTZ) was systematically lower 
than the energies of three plausible stack configurations. In order to get a better 
description of the van der Waals interaction, the density functional theory (DFT) 
results were adjusted by adding a phenomenological Lennard-Jones-type term into 
the total energy of the system. However, the van der Waals correction is some- 
what arbitrary and its magnitude can not be rigorously justified. Depending on the 
choice of the parameters in the phenomenological term both the T-shaped and the 
parallel-displaced (PD) stack configurations can be the global energy minimum of 
the system. A simple model is proposed which is useful for qualitative understanding 
of possible geometries of the coronene dimer and larger coronene clusters. The model 
represents coronene dimer as two sets of charged rings interacting via Coulomb and 
Lennard- Jones potentials. The model provides an intuitively clear explanation why 
the T-shaped dimers can be of importance even for some of moderately large PAHs 
such as coronene and, may be, for circumcoronene. The unexpectedly strong bind- 
ing of the coronene dimer in the T-shaped configuration is connected to the round 
shape of the coronene molecules. Indeed, rotation or parallel-displacement of the 
non-round monomers results in significantly smaller Coulomb repulsion as compared 
to the "face to face" sandwich configuration. On the other hand, rotation and/or 
displacement are much less effective for the stacked round molecules. Therefore, the 
round shape of the coronenes leads to an increased role of the electrostatic repulsion 
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in the stack configurations and the T-shaped configuration can become competi- 
tive. The proposed model can be easily generalized to other polycyclic aromatic 
hydrocarbons, DNA bases, etc. 

PACS numbers: 



I. INTRODUCTION 

Molecules of polycyclic aromatic hydrocarbons and their ions are of current interest as 
they are believed to be one of the most probable species in the interstellar space responsible 



stars into the intensive IR radiation, known as unidentified 
. The physical mechanisms and effectiveness of energy 



for converting the UV radiation o: 
infrared bands 

transfer from the UV to the IR radiation occuring via transient heating j8[ are currently 

being studied, see, e.g., liiy,Hiii|. 

Astrophysical observations suggest that the interstellar medium contains a variety of small 
carbonaceous part.cles wth the s.es ranging from dust grains down to singlemoleeules Q. 
The most abundant particles contain between 100 and 1000 carbon atoms [IJ] in the forms 
of collections (clusters) of PAH-like molecules [l^ of various degrees of hydrogenation P| . 
Hence, studying properties of PAH clusters presents a necessary step in further research of 
interstellar medium. 

The current knowledge of structure and ways of formation of PAH clusters is very limited, 
jaboratory studies have just started to appear and only anthracene Q, Q], naphthalene 
14 1 and coronene [3, ll^ clusters have been produced so far, so that the array of available 
experimental data is sparse. Theoretical investigations of oligomers of the larger PAHs are 
also scarce and, furthermore, their results are very sensitive to a subtle interplay between 
the van der Waals polarization forces (also sometimes called dispersive interaction) and 



the electrostatic interaction. Even for benzene, the simplest aromatic hydrocarbon, an 



*Electronic address: o.obolensky@fias.uni- frankfurt.de . On leave from: A.F. loffe Institute, Politechnich- 
eskaja str. 26, St. Petersburg f94021, Russia 

^^Electronic address: solovyov@fias.uni-frankfurt.de, On leave from: A.F. loffe Institute, Politechnicheskaja 
str. 26, St. Petersburg 194021, Russia 



3 



unequivocal choice between the T-shaped and the parallel-displaced stack configurations 
can not be made since the calculated binding energies for these configurations strongly 
depend on the choice of the method of calculation and on the basis set and they do not 
seem to be converging even at the most advanced levels of theory . Alternative 



approaches based on mixing of the ab initio and the phenomenological treatments 



lave been 
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mm), 



suggested which allow a more fiexible tuning of interactions (see, e.g., 
but their universality and predictive power are not quite clear yet. 

We have performed density functional theory calculations of energies for four plausible 
configurations [2^ of the coronene dimer. These included three stack configurations (super- 
imposed stack - "sandwich" (SS), twisted stack (TS) and parallel-displaced (PD) stack) and 
the T-shaped configuration. To our knowledge, this is the first comparative study of the 
configurations of the coronine dimer within the ab initio approach. The only other ab initio 
calculation of coronene dimer parameters was done for the PD stack configuration |2l| . 

The energy of the T-shaped configuration at all considered levels of theory (combinations 
of the B3LYP and PBE functionals with 6-31+G(d), D95, cc-pVDZ, and cc-pVTZ basis 
sets) was systematically lower than the energies of the stack configurations. In order to get 
a better description of the van der Waals interaction, the DFT results were adjusted by 
adding a phenomenological Lennard- Jones-type term into the total energy of the system. 
The problem of this approach is that the van der Waals correction is somewhat arbitrary and 
different authors choose different parameters even for the same systems (cf., e.g., j2l| and 
[2^). Depending on magnitude of the van der Waals term, both the T-shaped and the stack 
configurations can be the global minimum configuration of the system. Even though the 
definitive conclusion can not be drawn at the moment, it is quite surprising on its own that 
the T-shaped configuration is energetically competitive for such a large system for which a 
graphite-like stacked structure was expected [2l|. 

We propose a simple model useful for quahtative understanding of dimerization and clus- 
terization of polycyclic aromatic hydrocarbons. The model represents a PAH molecule as 
a set of circular or deformed rings corresponding to the chains of carbon atoms surrounded 
by a chain of hydrogens. For example, the coronene molecule is represented as a set of four 
nested circular rings with the outer ring corresponding to the hydrogens. The electron den- 
sity in PAHs is shifted from the hydrogens towards the interior of the molecules. In a crude 
approximation this charge distribution can be described as a polarized band (with a posi- 
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tively charged outer edge and a negatively charged inner edge) located along the perimeter 
of a molecule. This charge distribution can also be considered as a dipole distributed along 
the border of a molecule. Two PAH molecules interact electrostatically and via the van der 
Waals forces between the sets of the rings substituting the monomers. 

The model provides a qualitative explanation why even for such a large PAH as coronene 
(C24H12) the T-shaped configuration is comparable in energy with the stack configuration. 
This happens due to the round shape of the coronene molecule which leads to an increased 
electrostatic repulsion in the stack configurations. For non-round molecules in the stack 
configurations, the electrostatic interaction can be made less repulsive (or even attractive) 
by rotating or by parallel displacing the molecules in the dimer, while for stacked round 
molecules the rotation or displacement are much less effective and the T-shaped configuration 
can become competitive. 



II. DENSITY FUNCTIONAL THEORY TREATMENT OF CORONENE DIMER 

We have performed density functional theory calculations of binding energies for four 
configurations of coronene dimer. The dimer binding energy is defined as follows: 



2£'jnono -^dimer; (l) 

where -Edimer is the total energy of the dimer and -Emono is the energy of the single molecule. 
Note that according to this definition the binding energies of the stable bound states of the 
dimer are positive. 

The studied geometries are shown in Figure^ and the binding energies are summarized in 
TablelU The calculations have been done with the B3LYP and PBE functionals by expanding 
the molecular orbitals into the standard 6-31-|-G(d), D95, cc-pVDZ, and cc-pVTZ basis sets 

nn 

as implemented in Gaussian03 package |2fil. 12 7|. 

The energy of the T-shaped configuration at all considered levels of theory was system- 
atically lower than the energies of the three stack configurations. It is a quite surprising 
result for such a large PAH molecule. For the smallest aromatic hydrocarbon, benzene, the 
T-shaped configuration is almost isoenergetic with the parallel-displaced stack configuration 
[iSlllSl^II- It is commonly believed that increase in number of aromatic rings favors stack 

nnn 

configurations over the T-shaped ones |22, |2J, |2SI|. In accord with this assumption, the 
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TABLE I: The binding energies of the coronene dimer (in kcal/mol) calculated at various levels 
of theory. For notations and geometries of the dimer configurations see Figure ^ 





SS 


TS 


PD 


T 


B3LYP/6-31+G(d) 


-12.55 


-9.16 


-9.41 


0.33 


B3LYP/D95 


-2.99 


-2.64 


-1.91 


0.43 


B3LYP/cc-pVDZ 


-4.49 


-2.72 


-2.55 


0.43 


PBE/6-31+G(d) 


-5.34 


-0.91 


1.7 


1.98 


PBE/D95 


0.88 


1.39 


2.32 


2.83 


PBE/cc-pVDZ 


0.12 


0.66 


1.72 


1.79 


PBE/cc-pVTZ 


-0.51 


0.76 


1.22 


1.48 


PBE/TZV(2d,2p) from Ref. \21] 






-4.89 




BLYP/TZV(2d,2p) from Ref. [21j 






-17.52 





T-shaped (or, more generally, "edge to plane") dimers of pyridine, naphtalene, azulene and 
anthracene have been found less binding than the stacked c onfig urations j3, 24 , 13, ■ 
In the only other ab initio calculation of a coronene dimer j2l| a graphite-like (parallel- 
displaced) stacked structure was considered only. 

To our knowledge, this paper is the first comparative study of the possible configurations 
of the coronine dimer performed within the ab initio approach. There were also studies 
of the —e dimer which used phe„„,„e„o.og,ca. approaches WU. In the most 
recent study of this type j23|, the dimers and clusters of PAH molecules were modeled 
within a rigid-body approximation with the use of explicit interatomic potentials (a sum of 
an electrostatic and a van der Waals terms). None of these semiempirical studies found the 
T-shaped form of the coronene dimer to be energetically preferable. 

It is well known that the pure DFT methods does not reproduce the attractive polar- 
ization interactions well enough, for a discussion see, e.g., Ref. [2^ and references therein. 
This may present a serious problem in such intricate cases as determining the binding en- 
ergy for loosely bound dimers of PAH molecules. In order to estimate the influence of the 
polarization van der Waals interaction on the preferred dimer geometry, we have included a 
phenomenological Lennard- Jones-type term into the total energy of the system, similar to 




PD Y r T 





1^ 



FIG. 1: (Color online) Geometries of four configurations of the coronene dimer. In the sandwich 
(SS) configuration the distance R between the centers of mass of the monomers equals 4.27 A. In 
the twisted sandwich (TS) configuration R = 4.21 A and the rotation angle is 30°. In the parallel- 
displaced (PD) configuration R = 3.95 A and the shift is 2.70 A. In the T-shaped (T) configuration 
R = 8.02 A. The binding energies, calculated at various levels of theory, can be found in Table 



the procedure developed in 



2l| 



E = Edft + II -# - 



(2) 



Here -Edft is the total energy of the system calculated within the frame of the DFT. Index 
i enumerates atoms in one monomer, index j enumerates atoms in the other monomer; 
correspondingly, rij is the distance between atom "i" in the first monomer and atom "j" in 
;he second monomer. The coefficients Cg-' and C/2 can be tabulated directly (as is in Ref. 
3^ ) or found with the help of a combining rule (as is in Ref. 21 1) 



M 



M = 6, 12. 



(3) 



Unfortunately, there are no fixed values for the coefficients Cg and C12 which would be 
universally applicable in a wide scope of situations. Even for the same systems different 
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authors choose different parameters (cf., e.g., 21| and 2M)- The arbitrariness of the van 
der Waals correction seems to be a general problem of all such hybrid approaches. 

In order to illustrate how sensitive the overall conclusions may be with respect to the val- 
ues of the phenomeno logical parameters we present in Table ITTl the binding energies for each 
of the considered dimer configurations for various reasonable choices of the parameters Cq, 
C|^, Ci2, and C^. Depending on magnitude of the van der Waals term, both the T-shaped 
and the stack configurations can be the global minimum configuration of the system. Thus, 
we suppose that a definitive conclusion can not be drawn at the moment. An experimental 
measurement of certain characteristics (such as frequencies of transitions between rotational 
levels) would be useful for determining the global energy minimum configuration of the 
dimer. To facilitate identification of the dimer configurations, in Table IIIII the rotational 
constants for each configuration are given. The rotational constants -B.^, connected to the 
principal values of the inertia tensor 1 as Bj^ = h/ (47rX^^) (77 = x, y, z), define the rotational 
levels En = B^Ji^J + 1) corresponding to different rotational quantum numbers J. 



TABLE II: The binding energies (in kcal/mol) for four configurations of the coronene dimer (see 
Figure ^ calculated according to expressions ^ and for several sets of the coefficients Cq 
(in J/mol nm^) and C12 (in mJ/mol nm^^). The energies of the van der Waals terms are given 
in parentheses. The optimization of the dimer geometries and calculations of the DFT binding 
energies have been done with the PBE/D95 method. 





SS 


TS 


PD 


T 


= 2.77 Cf=O.U Cf-H = 0.63" 
= 4.87 = 0-09 (^12^ = 0-68 


25.73 (24.85) 


19.07 (17.68) 


17.72 (15.40) 


7.54 (4.71) 


C2 = 1.15 Cf = 0.11 ' 


17.97 (17.09) 


10.12 (8.73) 


9.97 (7.65) 


5.52 (2.69) 


C2 = 0.88 Cf = 0.11 
Cfj = 4.87 (7^2 = 0.09= 


2.34 (1.46) 


5.73 (4.34) 


6.00 (3.68) 


3.51 (0.68) 


C12 = 


0.88 (0.00) 


1.39 (0.00) 


2.32 (0.00) 


2.83 (0.00) 



TromRef. S2|. 

^From Ref. 2l| (Ce is adjusted by the factor 0.7 ir jm one of the values (1.65) given in Ref. (33i| ) . 
"A dump function has been used instead of C12 ,21] 

''Another possible value (1.26) of the coefficient Cg from Ref. _33J multiplied by 0.7 as in Ref. 
'^Coefficients C12 are from Ref. [3 



We note that the superimposed stack (sandwich) configuration probably is not a minimum 
on the potential energy surface. An analysis done with the use of a phenomenological 
potential has identified it as a transitional state between the twisted stack and the parallel- 
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displaced stack configurations. This conclusion is indirectly supported by the results of the 
high-level ab initio calculations which demonstrated that in benzene the superimposed stack 
is also a saddle point rather than a true minimum . 

TABLE III: Rotational constants (in GHz) for different configurations of the coronene dimer cal- 
culated for the geometries optimized at B3LYP/D95 and PBE/D95 levels of theory. For notations 
and geometries of the dimer configurations see Figure ^ 





SS 


TS 


PD 


T 


B3LYP/D95 


0.085 


0.085 


0.082 


0.087 


0.087 


0.082 


0.123 


0.060 


0.048 


0.110 


0.041 


0.036 


PBE/D95 


0.087 


0.087 


0.082 


0.089 


0.089 


0.082 


0.113 


0.074 


0.059 


0.111 


0.041 


0.037 



III. RINGS MODEL OF POLYCYCLIC AROMATIC HYDROCARBONS 

The model approaches sometimes can be very useful in revealing the underlying physics 
of a process and in pointing to the directions of search. In some situations, when the ab 
initio calculations are not feasible, the model calculations are the only approach one has to 
resort to. We present here a model which elucidates the role of the shapes of PAH molecules 
in the process of dimerization and in formation of bigger PAH clusters. 

The model represents PAH molecules as a set of circular or deformed rings which corre- 
spond to the chains of carbon atoms surrounded by a chain of hydrogens. For example, the 
coronene molecule can be described as consisting of four concentric circular rings represent- 
ing three chains of carbon atoms and one chain of hydrogens, see Figure |21 

The electron density in PAHs is shifted from the hydrogens towards the interior of the 
molecules. Therefore, the outer ring (corresponding to the hydrogens) is always charged 
positively. The inner rings can, in principle, bear any charge, but their net charge, in the 
case of a neutral molecule, must be negative and compensate the positive charge of the 
outermost ring. 

The magnitude of charge for each ring is defined on the basis of the results of the DFT 
calculations as the sum of the partial charges of the atoms belonging to the corresponding 
chain. Various methods can be used for assigning the partial charges to the atoms. For 
example, Mulliken suggested to calculate partial charges as the integral over the electron 
density contained within a sphere of a given radius, specific to the corresponding type of 
atoms js^. Alternatively, one can expand the electron density into a set of "natural bond 
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FIG. 2: (Color online) The round-shaped coronene molecule represented as a set of four nested 
circular rings. The radii of the rings are 1.438, 2.876, 3.595, and 4.691 A (Rc-c=l-438A, 
Rc-H=l-096A). The outermost ring corresponds to the hydrogens and is positively charged (Q=1.8 
e). In the simplest case, this charge is compensated by the negative charge of the next ring, cor- 
responding to the chain of outer carbon atoms. Two inner rings, then, are neutral. 



orbi 

m 



;als" (NBO) chosen in such a way that each orbital is occupied by, ideally, two electrons 



36 



We use the so-called Electrostatic potential scheme (ESP), in which the partial 
charges are determined from the best fit of the actual electrostatic potential of the molecule 
j^]. This method is not only better justified physically for our purposes, but it is also less 
dependent on the choice of basis set used in the calculations, see Table IIVI 

In the simplest approximation, the electron density, shifted from the hydrogens, is as- 
signed to the outer carbons, covalently bonded to the hydrogens. Then the ring correspond- 
ing to the chain of the outermost carbons is charged negatively. The inner rings remain 
neutral. The charge distribution in such approximation can be crudely described as a po- 
larized band (with a positively charged outer edge and a negatively charged inner edge) 
located along the perimeter of the molecule. This charge distribution can also be considered 
as a dipole continuously distributed along the outer border of the molecule. In fact, within 
this approximation one can assign a typical value of the partial charge (say, 0.15 j32|) to all 
the hydrogen atoms and then one immediately obtains the charges of the two outer rings. 
The DFT analysis of the electron density in this case is not necessary. This approximation 
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TABLE IV: Partial charges in coronene assigned to carbon and hydrogen atoms accordin g to 
the Muhiken natural bond orbitals (NBO) [s^ and best fit electrostatic potential (ESP) |27| 
schemes (in units of elementary charge). Ci are carbons from the innermost chain (there are 6 
of them), C2 are carbons from the intermediate carbon chain (6), C3 are carbons from the outer 
carbon chain (12). 





MuUiken 


NBO 


ESP 




Ci 


C2 


C3 


H 


Ci 


C2 


C3 


H 


Ci 


C2 


C3 


H 


B3LYP/6-31+G{d) 


-0.346 


0.741 


-0.384 


0.185 


0.158 


-0.138 


-0.213 


0.241 


-0.002 


0.129 


-0.207 


0.148 


B3LYP/D95 


-0.008 


0.463 


-0.46 


0.232 


-0.016 


-0.041 


-0.19 


0.218 


0.002 


0.13 


-0.237 


0.168 


B3LYP/cc-pVDZ 


0.064 


-0.02 


0.026 


-0.049 


-0.008 


-0.053 


-0.199 


0.229 


-0.004 


0.095 


-0.171 


0.129 


PBE/6-31+G(d) 


-0.362 


0.807 


-0.42 


0.197 


0.115 


-0.136 


-0.22 


0.251 


-0.003 


0.121 


-0.223 


0.148 


PBE/D95 


0.033 


0.454 


-0.482 


0.238 


-0.016 


-0.044 


-0.196 


0.225 


-0.003 


0.133 


-0.236 


0.16 


PBE/cc-pVDZ 


0.076 


-0.014 


0.009 


-0.04 


-0.008 


-0.055 


0.207 


0.239 


-0.0002 


0.086 


-0.18 


0.13 


PBE/cc-pVTZ 


0.371 


-0.029 


-0.335 


0.164 


-0.009 


-0.054 


-0.187 


0.218 


0.002 


0.114 


-0.205 


0.149 



provides a tool for studying interactions of large PAHs which cannot be treated with ab 
initio methods. 

The interaction energy between two PAH molecules consists of a Coulomb and a van der 
Waals terms calculated for each pair of the rings: 



(4) 



a,l3 



where index a enumerates the rings in one monomer, and index /3 enumerates the rings in 
the other monomer. 

The energy of Coulomb interaction of two charged rings i?^"' can be found by integrating 
the electrostatic potential created by one ring along the contour of the other ring: 



Pi 



(5) 



where QpSs the charge of the ring (3 equal to the sum of partial charges of all atoms belonging 
to the ring, is the length of the ring. The electrostatic potential 0^°"'(r) created by a 
circular ring is expresses via the complete Elliptic function 



r) 



2Q„ 1 1 



K 



TT ^2rRa_ + 1 Vx + l 



X 



2rRn sin 9 



(6) 
(7) 
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Here is the radius of the ring a, r is the modulus of r, 6 is the angle between r and the 
ring normal. 

The van der Waal's term E^^^ can be calculated similarly, 

- f (g^r(r) + g^J(r)) M, (8) 

The van der Waals charges and are defined as sums of the Cg and C12 coefficients, 
respectively, of all the atoms belonging to the ring. For a circular ring the van der Waals 
potentials 0^ and (fi^ are expressed via elementary functions. 



Qi^ 1 -I- 1 



2 (2ri?„)3 (x2 - 1)5/2^ 



X. . Ql 1 x(8x^ + 40x^ + 15) 

^a(r) = ^ /o^D .a 7-0 TTTTT^ • (10) 



8 {2rR^Y 1)11/2 
The model predicts the T-shaped and PD configurations to be stable, in accordance with 
the results of the DFT calculations. In addition to it, the model provides an intutively 
obvious explanation why even for such a large PAH as coronene the T-shaped configura- 
tion is comparable in energy with the stack configuration. This happens due to the round 
shape of the coronene molecule, since only the "edge to plane" configurations can be stable 
electrostatically. On the other hand, the dispersive interaction favors superimposed stacks 
for which the "contact area" is largest. For the non-round molecules, e.g., anthracene, the 
interplay between these two factors results in parallel-displaced or rotated configurations 



24| , because a parallel displacement or a rotation can significantly decrease the electrostatic 



repulsion while not seriosly affecting the van der Waals attraction. However, for the round 
molecules, such as coronene, the rotation or parallel displacement are much less effective in 
decreasing the Coulomb repulsion and the T-shaped configuration can become competitive. 

Therefore, we conclude that in the dimerization process (and in formation of larger clus- 
ters of PAH molecules) the shapes of the molecules are very important. The proposed simple 
model can be useful in qualitative analysis of the possible configurations and in making the 
results of such analysis intuitively clear. The model can be applied to any polycyclic aro- 
matic hydrocarbons, and it can be easily generalized on the DNA bases, etc. The model can 
also be used for predicting structures of larger coronene clusters and clusters of other PAHs. 
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